For generating most, but not all figures in the manuscript.
library(tidyverse)
library(cowplot)
library(lubridate)
library(mgcv)
source("UVP_2017_library.R")
Parsed with column specification:
cols(
Cruise = [31mcol_character()[39m,
Station = [31mcol_character()[39m,
`mon/dd/yyyy` = [31mcol_character()[39m,
`hh:mm` = [34mcol_time(format = "")[39m,
`Longitude [degrees east]` = [32mcol_double()[39m,
`Latitude [degrees north]` = [32mcol_double()[39m,
`Bottom Depth [m]` = [32mcol_double()[39m,
`Pressure [db]` = [32mcol_double()[39m,
`Temperature [degrees C]` = [32mcol_double()[39m,
`Temperature 2 [degrees C]` = [32mcol_double()[39m,
`Salinity [psu]` = [32mcol_double()[39m,
`Salinity 2 [psu]` = [32mcol_double()[39m,
`Fluorescense [mg/m^3]` = [32mcol_double()[39m,
`Beam Transmission [%]` = [32mcol_double()[39m,
PAR = [32mcol_double()[39m,
`Oxygen [umol/kg]` = [32mcol_double()[39m,
`Oxygen [% saturation]` = [32mcol_double()[39m
)
theme_set(theme_cowplot())
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
These read back in in UTC, since write_csv saves everthing into that time zone.
options(readr.default_locale=readr::locale(tz="Mexico/General"))
bes<- read_csv("dataOut/binned_EachSize.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
bds <- read_csv("dataOut/binned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
ues <- read_csv("dataOut/unbinned_EachSize.csv")
Parsed with column specification:
cols(
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m,
depth = [32mcol_double()[39m,
psd_gam = [32mcol_double()[39m,
vol = [32mcol_double()[39m,
sizeclass = [31mcol_character()[39m,
lb = [32mcol_double()[39m,
ub = [32mcol_double()[39m,
binsize = [32mcol_double()[39m,
TotalParticles = [32mcol_double()[39m,
nparticles = [32mcol_double()[39m,
n_nparticles = [32mcol_double()[39m,
biovolume = [32mcol_double()[39m,
speed = [32mcol_double()[39m,
flux = [32mcol_double()[39m,
flux_fit = [32mcol_double()[39m,
GamPredictTP = [32mcol_double()[39m
)
uds <- read_csv("dataOut/unbinned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
Specify the base of the photic zone (which is inside of the OMZ) and the bae of the OMZ
PhoticBase <- 160
OMZBase <- 900
DVMBase <- 600
OMZTop <- 90
library(scales)
#https://stackoverflow.com/questions/30179442/plotting-minor-breaks-on-a-log-scale-with-ggplot
log_breaks = function(maj, radix=10) {
function(x) {
minx = floor(min(logb(x,radix), na.rm=T)) - 1
maxx = ceiling(max(logb(x,radix), na.rm=T)) + 1
n_major = maxx - minx + 1
major_breaks = seq(minx, maxx, by=1)
if (maj) {
breaks = major_breaks
} else {
steps = logb(1:(radix-1),radix)
breaks = rep(steps, times=n_major) +
rep(major_breaks, each=radix-1)
}
radix^breaks
}
}
scale_x_log_eng = function(..., radix=10) {
scale_x_continuous(...,
trans=log_trans(radix),
breaks=log_breaks(TRUE, radix),
minor_breaks=log_breaks(FALSE, radix))
}
ybreaks = seq(from = 1200, to = 0, by = -50)
ylabels = ybreaks
ylabels[c(FALSE, TRUE, TRUE, TRUE)] <- ""
#theme_set(theme_bw)
#theme_set(theme_cowplot)
PlotParticlesmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = tot_nparticles, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse(limits = c(1200, 0),
breaks = ybreaks,
labels = ylabels) +
scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black"), limits = c(0, 24)) +
scale_x_log10(breaks = c(10, 100, 1000), minor = c(5, 50, 500)) +
#theme(legend.position = "none") +
#scale_x_log_eng()+
labs(y = "Depth (m)", x = "Particles / L") +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue") +
geom_hline(yintercept = OMZTop, color = "darkblue") +
#geom_hline(yintercept = DVMBase, color = "darkgoldenrod", lty = "dashed") +
theme(legend.position = "none")
PlotPSDmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = psd, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse(limits = c(1200, 0), breaks = ybreaks, labels = ylabels) +
scale_shape_manual(name = "Day", values = c(21:25)) +
scale_fill_gradientn(name = "Hour",breaks = c(0, 6, 12, 18, 24), labels = c("0", "6", "12", "18", "24"), colors = c("black", "blue", "white", "orange", "black"), limits = c(0, 24)) +
labs(y = "Depth (m)", x = "Particle Size Distribution Slope") +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue") +
geom_hline(yintercept = OMZTop, color = "darkblue") +
#geom_hline(yintercept = DVMBase, color = "darkgoldenrod", lty = "dashed") +
coord_cartesian(clip = "off") +
annotate(geom = "text", x = -2, y = 1100, label = "More large \n particles", inherit.aes = FALSE, size = 6) +
annotate(geom = "text", x = -4.5, y = 1100, label = "More small \n particles", inherit.aes = FALSE, size = 6)
Ignoring unknown parameters: inherit.aesIgnoring unknown parameters: inherit.aes
# annotation_custom(
# grob = textGrob(label = "More big particles", hjust = 0, gp = gpar(cex = 12)),
# ymin = 0, # Vertical position of the textGrob
# ymax = 0,
# xmin = 0, # Note: The grobs are positioned outside the plot area
# xmax = 0)
plot_grid(
PlotParticlesmany,
PlotPSDmany,
rel_widths = c(1, 2),
labels = c("A", "B")
)
Removed 266 rows containing missing values (geom_point).Removed 266 rows containing missing values (geom_point).
ggsave("figures/ParticlesPSDMany.png")
Saving 12 x 7.41 in image
ggsave("figures/ParticlesPSDMany.svg")
Saving 12 x 7.41 in image
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(tot_nparticles~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG2 <- gam(tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG3 <- gam(tot_nparticles ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
#FSG4 <- gam(tot_nparticles~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4,
bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.96538 0.09655 92.86 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.0000 1.000 5.663 0.0202 *
s(Day) 1.4473 1.694 2.033 0.0921 .
s(Hour) 0.6402 2.000 0.502 0.2154
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.135 Deviance explained = 17.4%
GCV = 0.68369 Scale est. = 0.64319 n = 69
summary(FSG2)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.96538 0.09727 92.17 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.000 1.000 5.741 0.0194 *
s(Day) 1.469 1.718 2.214 0.0792 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.122 Deviance explained = 15.3%
GCV = 0.68744 Scale est. = 0.65288 n = 69
summary(FSG3)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.9654 0.1004 89.34 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1 1 5.736 0.0194 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.0651 Deviance explained = 7.89%
GCV = 0.71559 Scale est. = 0.69485 n = 69
#summary(FSG4)
summary(FSG1)$r.sq - summary(FSG2)$r.sq # extra R^2 explained by hour
[1] 0.01304448
summary(FSG2)$r.sq - summary(FSG3)$r.sq # extra explained by day
[1] 0.05646094
summary(FSG3)$r.sq # R^2 explained by depth
[1] 0.06511156
There is no statisticlly significant affect of time on particle number. If I take all of the time variables out and just compare to depth, there is a relationship to depth p = 0.02. But the R^2 is only 6.5%. Pretty weak
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(psd~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG2 <- gam(psd ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG3 <- gam(psd ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
FSG4 <- gam(psd~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.90483 0.01809 -215.8 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.848 1.977 88.789 <2e-16 ***
s(Day) 1.506 1.756 0.953 0.4853
s(Hour) 1.361 2.000 2.260 0.0424 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.716 Deviance explained = 73.6%
GCV = 0.024631 Scale est. = 0.022591 n = 69
#summary(FSG2)
#summary(FSG3)
summary(FSG4)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.90483 0.01814 -215.2 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.845 1.976 88.630 <2e-16 ***
s(Hour) 1.524 2.000 2.559 0.0409 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.714 Deviance explained = 72.8%
GCV = 0.024244 Scale est. = 0.022709 n = 69
summary(FSG3)$r.sq # R^2 of depth
[1] 0.6929817
summary(FSG4)$r.sq - summary(FSG3)$r.sq # Improvement from adding hour of the day
[1] 0.02133449
summary(FSG1)$r.sq - summary(FSG4)$r.sq # Improvement from then adding day of the week
[1] 0.001481684
PSD varies with depth, but doesn’t statistically relate to hor orday. Comparing the R2 values from models tells us that you explain 69% of varience with depth.
PlotNParticlesEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = tot_nparticles, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
#geom_path(aes(x = tot_nparticles)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + scale_color_manual(values = c("gray20", "brown")) +
labs(x = "Particles/L", y = "Depth (m)") +
theme(legend.position = "none") +
scale_shape_manual(values = c(1:5)) +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = 200, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue")
PlotNParticlesEP
I removed one outlyer from p16 for visualization purposes (300 particles/l at surface)
PlotPSDEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = psd, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
geom_path(aes(x = psd_gam)) +
geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_color_manual(values = c("gray20", "brown")) +
scale_shape_manual(values = c(1:5)) + labs(y = "", x = "Particle Size Distribution Slope") +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = 200, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue")
PlotPSDEP
plot_grid(PlotNParticlesEP, PlotPSDEP, rel_widths = c(2,3), labels = c("A", "B"))
Removed 1211 rows containing missing values (geom_point).Removed 1211 rows containing missing values (geom_point).Removed 1211 row(s) containing missing values (geom_path).
ggsave("figures/ParticlesAndPSD_ETNPVsP16.svg")
Saving 10 x 4 in image
ggsave("figures/ParticlesAndPSD_ETNPVsP16.png")
Saving 10 x 4 in image
Large and spall particle number, flux and size
mainParticleComponents <- bds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
select(project, profile, depth,
tot_nparticles, small_nparticles, big_nparticles,
tot_psd = psd, small_psd, big_psd,
tot_flux_fit, small_flux_fit, big_flux_fit) %>%
pivot_longer(cols = -c("project", "profile", "depth")) %>%
separate(name, c("size", "meas")) %>%
mutate(meas = recode(meas, nparticles = "particles/L")) %>%
mutate(meas = factor(meas, levels = c("particles/L", "flux", "psd")))
Expected 2 pieces. Additional pieces discarded in 273 rows [7, 8, 9, 16, 17, 18, 25, 26, 27, 34, 35, 36, 43, 44, 45, 52, 53, 54, 61, 62, ...].
PlotFlx <- mainParticleComponents %>%
filter(meas != "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size ~ meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + theme(axis.title.x = element_blank(), legend.position = "none", strip.background.y = element_blank(), strip.text.y = element_blank(), plot.margin = unit(c(7,0,7,7), "pt")) + scale_color_manual(values = c(P16 = "brown", ETNP = "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = PhoticBase, color = "darkgreen")
PlotPSD <- mainParticleComponents %>%
filter(meas == "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size~meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.line.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), plot.margin = unit(c(7,7,26.5,0), "pt")) +
scale_color_manual(values = c(P16 = "brown", ETNP = "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = PhoticBase, color = "darkgreen")
plot_grid(PlotFlx, PlotPSD, rel_widths = c(3, 2))
Removed 246 rows containing missing values (geom_point).Removed 123 rows containing missing values (geom_point).
ggsave("figures/BigVsSmall.svg")
Saving 7.29 x 4.5 in image
ggsave("figures/BigVsSmall.png")
Saving 7.29 x 4.5 in image
Flux small and flux tot track so closely because particle fractal dimension alpha, plus flux fractal dimension, gamma > |psd|. since the size distribution of the flux sould be PSD + ag (psd is negative in this case). Yo ucan see the variance at the one depth where psd is flatest at the very top.
Example particle size distributions
eg_dataline <- bds %>%
filter(profile == "stn_043", depth == 162.5)
eg_slope = eg_dataline %>% pull(psd)
eg_icp = eg_dataline %>% pull(icp)
eg_vol = eg_dataline %>% pull(vol)
eg_datablock <- bes %>%
filter(profile == "stn_043", depth == 162.5)
eg_lb = eg_datablock$lb
eg_binsize = eg_datablock$binsize
eg_nnp = exp(eg_icp + log(eg_lb) * eg_slope)
eg_np = eg_nnp * eg_binsize
eg_tp = eg_np * eg_vol
eg_df <- tibble(lb = eg_lb, n_nparticles = eg_nnp, nparticles = eg_np, TotalParticles = eg_tp)
EgNNP <- eg_datablock %>%
ggplot(aes(x = lb, y = n_nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Binsize & Volume Normalized \n Particles (#/L/mm)", x = "Size (mm)")
EgNP <- eg_datablock %>%
ggplot(aes(x = lb, y = nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Normalized Particles" , x = "Size (mm)")
EgTP <- eg_datablock %>%
ggplot(aes(x = lb, y = TotalParticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs( y = "Total Particles Observed (#)", x = "Size (mm)")
plot_grid(EgNNP, EgTP, labels = c("A", "B"))
Transformation introduced infinite values in continuous y-axisTransformation introduced infinite values in continuous y-axis
ggsave("figures/ExamplePSD163m.png")
Saving 7.29 x 4.5 in image
ggsave("figures/ExamplePSD163m.svg")
Saving 7.29 x 4.5 in image
Flux attenuation with respect ot depth and time. All extrapolated from the UVP and traps combined.
scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_b <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_c <- function(x) {
xout <- gsub("1e", "10^{", format(x),fixed=TRUE)
xout <- gsub("{-0", "{-", xout,fixed=TRUE)
xout <- gsub("{+", "{", xout,fixed=TRUE)
xout <- gsub("{0", "{", xout,fixed=TRUE)
xout <- paste(xout,"}",sep="")
return(parse(text=xout))
}
scale_x_log10nice <- function(name=NULL,omag=seq(-10,20),...) {
breaks10 <- 10^omag
scale_x_log10(breaks=breaks10,labels=scientific_10_c(breaks10),...)
}
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlx <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(1000, 0), breaks = ybreaks, labels = ylabels) +
scale_x_log10nice()+
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black"), limits = c(0, 24)) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
#labs(x = "moo", y = "Depth (m)") +
geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 20, xmax = 180, ymin = 75, ymax = 500), colour = "red", fill = NA, inherit.aes = FALSE) +
theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm")) +
geom_segment(aes(y = PhoticBase, yend = PhoticBase, x = 20, xend = 500), color = "darkgreen", stroke = 0.5)+
geom_segment(aes(y = OMZBase, yend = OMZBase, x = 20, xend = 500), color = "darkblue", stroke = 0.5) +
geom_segment(aes(y = OMZTop, yend = OMZTop, x = 20, xend = 500), color = "darkblue", stroke = 0.5)
Ignoring unknown parameters: strokeIgnoring unknown parameters: strokeIgnoring unknown parameters: stroke
#geom_segment(aes(y = DVMBase, yend = DVMBase, x = 20, xend = 500), color = "darkgoldenrod", lty = "dashed", stroke = 0.5)
#+ geom_hline(yintercept = OMZBase, color = "darkblue")
pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
pltFlxLegend <- get_legend(pltFlx)
Removed 14 rows containing missing values (geom_point).
pltFlx
#plotly::ggplotly(plt1)
Zooming in on where the action is happening
ylabels_fine <- ybreaks
ylabels_fine[c(FALSE, TRUE)] <- ""
#xbreaks_fine <- c(seq(from = 20, to = 50, by = 10), seq(from = 60, to = 180, by = 20))
xbreaks_fine <- seq(from = 20, to = 180, by = 10)
xlabels_fine <- xbreaks_fine
toss <- c(rep(FALSE, 5), rep(c(TRUE, FALSE), 6))
xlabels_fine[toss] <- ""
xlabels_fine[xlabels_fine %in% c(120, 160)] <- ""
#lenx <- length(xlabels_fine)
#xlabels_fine[c(lenx, lenx-2)] <- ""
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxZoom <- bds %>% filter(project == "ETNP" & depth <= 500 & depth >= 75) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(500, 75), breaks = ybreaks, labels = ylabels_fine) +
#scale_x_log10() +
scale_x_log10(breaks = xbreaks_fine,
labels = xlabels_fine,
limits = c(20, 180)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(values = rep(21:25, 2)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Smoothed Flux", y = "Depth") + theme(legend.position = "none")+
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = OMZTop, color = "darkblue") +
theme(axis.text.x = element_text(vjust = .5))
pltFlxZoom
#plotly::ggplotly(plt1)
Rate of change of flux, taken to the fifth root so one can see patterns.
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3 <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(1000, 0), breaks = ybreaks, labels = ylabels) +
scale_x_continuous(limits = c(-2.1, .6), breaks = seq(from = -2, to = .75, by = 0.5)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")+
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue") +
geom_hline(yintercept = OMZTop, color = "darkblue")
#geom_hline(yintercept = DVMBase, color = "darkgoldenrod", lty = "dashed")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3
#plotly::ggplotly(plt1pos)
Within panel drawing
pgTop <- ggdraw(pltFlxNoLegend
) +
draw_plot(pltFlxZoom, .4, .25, .55, .60) +
draw_plot_label(
c("","B"),
c(.05, 0.55),
c(1, 0.85),
size = 16
)
Removed 14 rows containing missing values (geom_point).
pgTop
pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c(“C”, ""), label_size = 14)
I don’t know whats going on below here
pgBottom <- pltDelta3 + geom_rect(aes(xmin = -2, xmax = -1.15, ymin = 170, ymax = 1000), colour = "gray50", fill = "white", inherit.aes = FALSE) + draw_plot(pltFlxLegend , -1.9, -575, .7)
pgBoth <- plot_grid(pgTop + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
pgBottom + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
ncol = 1, rel_heights = c(4, 4), labels = c("A", "C"), label_size = 16)
Removed 33 rows containing missing values (geom_point).
pgBoth
ggsave("figures/FluxDeepDive.png")
Saving 5 x 9 in image
ggsave("figures/FluxDeepDive.svg")
Saving 5 x 9 in image
Test for day to day and hourly variability in rate of change of flux (fifth root transformed)
There is variability with respect to depth, and day and hour Depth p = 0.03 R^2 = 0.088. Add affect of day p = 0.004, extra R^2 = 0.11, Add affect of hour p = 0.02 extra R2 = 0.12
bdsAddTime <- bds %>%
mutate(Hour = hour(time) + minute(time)/60, Day = day(time) + hour(time)/24 + minute(time)/24/60)
DFG1 <- gam(pracma::nthroot(DF/DZ, 5)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG2 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG3 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG_DayOnly <- gam(pracma::nthroot(DF/DZ, 5) ~ s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG_NoFifth <- gam(pracma::nthroot(DF/DZ, 1)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
summary(DFG1)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3) +
s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.05998 -3.065 0.00407 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.672 1.893 3.955 0.0606 .
s(Day) 1.907 1.991 4.661 0.0190 *
s(Hour) 0.774 2.000 0.658 0.1961
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.254 Deviance explained = 33.4%
GCV = 0.1732 Scale est. = 0.15112 n = 42
summary(DFG_DayOnly)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.06485 -2.835 0.00722 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Day) 1.873 1.984 3.391 0.0404 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.129 Deviance explained = 16.8%
GCV = 0.18962 Scale est. = 0.17665 n = 42
summary(DFG_NoFifth)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 1) ~ s(depth, k = 3) + s(Day, k = 3) +
s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.028354 0.007924 -3.578 0.000979 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.809e+00 1.963 9.247 0.00138 **
s(Day) 1.827e+00 1.970 3.361 0.06222 .
s(Hour) 1.396e-08 2.000 0.000 0.71927
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.341 Deviance explained = 39.9%
GCV = 0.0029647 Scale est. = 0.0026375 n = 42
# summary(DFG2)
# summary(DFG3)
#
# summary(DFG1)$r.sq - summary(DFG2)$r.sq
# summary(DFG2)$r.sq - summary(DFG3)$r.sq
# summary(DFG3)$r.sq
Time is now continuous day and hour
Plot of the gams above
#plot.new()
FluxGamPlot <- function(){
par(mfrow = c(2,2))
plot(DFG1)
abline(h = 0, col = "gray30", lwd = 2)
mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,1))
abline(h = 0, col = "gray30", lwd = 2)
mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,2))
abline(h = 0, col = "gray30", lwd = 2)
mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
}
FluxGamPlot()
png(filename = "./figures/FluxGamPlot.png", width = 10, height = 8, units = "in", res = 200)
FluxGamPlot()
dev.off()
png
2
(u mol C / m^3 / day)
disagFig <- bds %>% filter(project == "ETNP") %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2) +
scale_y_reverse(limits = c(1000, 0), breaks = ybreaks, labels = ylabels) +
scale_x_continuous(limits = c(-1, 1)) +
geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
#labs(x = bquote("Observed - Modeled Small Particle Flux"~(μmol/m^3/day)), y = "Depth (m)") +
labs(x = paste("Deviation from Model", expression((μmol/m^3/day)))) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
geom_hline(yintercept = OMZBase, color = "darkblue") +
geom_hline(yintercept = OMZTop, color = "darkblue") #+
#geom_hline(yintercept = DVMBase, color = "darkgoldenrod", lty = "dashed")
disagFig
#ggsave("..figures/FluxSizeShift.svg"
ggsave("figures/FluxSizeShift.png")
Saving 6 x 4 in image
ggsave("figures/FluxSizeShift.svg")
Saving 6 x 4 in image
# bdsAddTime <- bds %>%
# mutate(Hour = hour(time), Day = day(time))
bdsAddTime <- bds %>%
mutate(Hour = hour(time) + minute(time)/60, Day = day(time) + hour(time)/24 + minute(time)/24/60)
OZG1 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
OZG2 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
OZG3 <- gam(ospsDZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))
summary(OZG1)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.08042 0.01049 7.665 1.33e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.584 1.827 19.185 4.77e-06 ***
s(Day) 1.771 1.946 8.926 0.00173 **
s(Hour) 1.266 2.000 1.964 0.05140 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.416 Deviance explained = 45.5%
GCV = 0.0082683 Scale est. = 0.0075948 n = 69
summary(OZG2)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.08042 0.01079 7.451 2.85e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.558 1.805 18.393 8.14e-06 ***
s(Day) 1.838 1.974 6.924 0.00415 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.382 Deviance explained = 41.3%
GCV = 0.0085853 Scale est. = 0.0080382 n = 69
summary(OZG3)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.08042 0.01169 6.877 2.59e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.453 1.7 17.31 2.23e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.274 Deviance explained = 29%
GCV = 0.0097831 Scale est. = 0.0094354 n = 69
summary(OZG1)$r.sq - summary(OZG2)$r.sq # Extra from Hour
[1] 0.03410746
summary(OZG2)$r.sq - summary(OZG3)$r.sq # Extrafrom Day
[1] 0.1074734
summary(OZG3)$r.sq # Depth
[1] 0.2741944
Plot of those gams Figure S10
OSMSGamPlot <- function(){
par(mfrow = c(2,2))
plot(OZG1)
par(mfg = c(1,1))
abline(h = 0, col = "gray30", lwd = 2)
mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,2))
abline(h = 0, col = "gray30", lwd = 2)
mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(2,1))
abline(h = 0, col = "gray30", lwd = 2)
mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
}
OSMSGamPlot()
png(filename = "./figures/OSMSGamPlot.png", width = 10, height = 8, units = "in", res = 200)
OSMSGamPlot()
dev.off()
png
2
trapFlux3 <- read_csv("dataOut/fluxMS_distilled.csv")
Parsed with column specification:
cols(
Class = [31mcol_character()[39m,
Depth = [32mcol_double()[39m,
TrapID = [31mcol_character()[39m,
TrapType = [31mcol_character()[39m,
SampleType = [31mcol_character()[39m,
C_flux = [32mcol_double()[39m,
C_flux_umol = [32mcol_double()[39m
)
UVPFluxComb <- read_csv("dataOut/CombinedProfileFluxEst_DS.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
Flux = [32mcol_double()[39m
)
UVPFluxOE <- read_csv("dataOut/ObservedVsExpectedFlux.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
tn_flux = [32mcol_double()[39m,
profile = [31mcol_character()[39m,
project = [31mcol_character()[39m,
time = [31mcol_character()[39m,
tot_flux2 = [32mcol_double()[39m
)
fluxMS_distilled_toPlot <- trapFlux3 %>%
mutate(SampleType = recode(SampleType, `plus.p` = "plus-particles", top = "top-collector"))
Remove traps where mass spec didn’t work correctly 2-17 150 1-12 73m 1-12 148 2-14 100 |(TrapID == “2-17” & Depth == 150)
fluxMS_distilled_toPlot2 <- fluxMS_distilled_toPlot %>%
filter(!((TrapID == "1-12") | (TrapID == "2-14" & Depth == 100)|(TrapID == "2-17" & Depth == 150)))
#fluxMS_distilled_toPlot2
Traps where mass spec didn’t work.
UVPFluxPlot00 <- UVPFluxComb %>%
ggplot(aes(y = depth)) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(0, 200)) +
geom_point(aes(y = Depth, x = C_flux_umol, shape = TrapType, ID = TrapID),
colour = "black", stroke = 1, size = 5, data = fluxMS_distilled_toPlot2) +
geom_line(aes(x = Flux), size = 1, color = "black") +
geom_point(aes(x = -1, y = -1, size = "UVP Estimate")) + # dummy point for the legend
geom_point(aes(x = tot_flux2), size = 3, shape = 21, color = "white", fill = "black", data = UVPFluxOE) +
scale_shape_manual(values = c(25, 22))+
scale_size_manual(values = 1, name = "") +
ylab("Depth (m)") +
#xlab(expression(Flux µmolC/m^2/day)) +
xlab(expression(paste("x axis ", ring(A)^2))) +
xlab(expression(paste("Flux (µ mol C/", m^2, "/day)"))) +
guides(fill = guide_legend(override.aes = list(shape = 21))) +
theme_cowplot() +
theme(
legend.position = c(0.5, 0.4),
legend.box.background = element_rect(color = "black", size = 0.5),
legend.margin = margin(-10, 5, 10, 5)
)
Ignoring unknown aesthetics: ID
# UVPFluxPlot <- UVPFluxPlot00 +
# geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 15, xmax = 32, ymin = 45, ymax = 195), colour = "red", fill = NA, inherit.aes = FALSE)
UVPFluxPlot00
ggsave("figures/FittedFlux.png")
Saving 7.29 x 4.5 in image
ggsave("figures/FittedFlux.svg")
Saving 7.29 x 4.5 in image
TPPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = TotalParticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "TotalParticles Observed (#)")
nnpPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = n_nparticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() +
#geom_vline(xintercept = 1) + geom_vline(xintercept = 5) +
labs(y = "Depth (m)", x = "Binsize and Volume Normalized Particles (#/L/mm)")
FitPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = nnp_smooth, xmin = nnp_lower, xmax = nnp_upper, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() +
#geom_vline(xintercept = 1) + geom_vline(xintercept = 5) +
labs(y = "Depth (m)", x = "Smoothed - Normalized Particles (#/L/mm)") + geom_errorbar(width = 10, alpha = 0.5)
npLegend <- get_legend(FitPlot + theme(legend.box.margin = margin(0, 0, 40, 200)) + labs(col = expression(log[e](Size (mm)))))
Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).
plot_grid(
TPPlot + theme(legend.position = "none"),
nnpPlot + theme(legend.position = "none"),
npLegend ,
FitPlot + theme(legend.position = "none"),
labels = c("A", "B", "", "C")
)
Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).
ggsave("figures/AllParticleSizes.svg")
Saving 10 x 6.18 in image
ggsave("figures/AllParticleSizes.png")
Saving 10 x 6.18 in image
We are smoothing the station 043 data, with respect to depth and time. We are using this station because it is the only one to extend past 2000m. We find that the other profiles seem to give strange values near the ends of the ranges of the gams.
SameGam <- gam(TotalParticles ~s(log(lb), log(depth)), offset = log(vol * binsize), family = nb(),
data = bes %>% filter(project == "ETNP", depth <= 2000, profile == "stn_043")) # Looks good!
gam.check(SameGam)
Method: REML Optimizer: outer newton
full convergence after 4 iterations.
Gradient range [1.180573e-06,5.441536e-06]
(score 2243.109 & scale 1).
Hessian positive definite, eigenvalue range [5.359006,135.6228].
Model rank = 30 / 30
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(log(lb),log(depth)) 29 21 0.61 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
besE <- bes %>% filter(project == "ETNP")
lb_new <- exp(seq(from = log(0.1), to = log(2.1), by = 0.05))
ub_new <- lead(lb_new)
binsize_new <- ub_new - lb_new
lbbs <- tibble(lb = lb_new, ub = ub_new, binsize = binsize_new)
Expanded <- expand_grid(lb = exp(seq(from = log(0.1), to = log(2), by = 0.05)), depth = seq(from = 20, to = 2000, by = 20), time = as.factor(unique(besE$time))) %>% left_join(lbbs, by = "lb")
Pred <- exp(predict(SameGam, Expanded))
ToPlot <- bind_cols(Expanded, nnparticles = Pred) %>% mutate(time = as.character(time)) %>% mutate(nparticles = nnparticles * binsize)
WBColorMap <- ToPlot %>% filter(lb <= 2) %>%
ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c(name = expression(log[10](Particles/m^3/mm))) + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = OMZBase, color = "darkblue") + labs(y = "Depth (m)", x = "Size (mm)")
WBColorMap
mbGam <- ToPlot %>% group_by(depth) %>% nest() %>%
mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>%
mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2]))
pWBPSD <- mbGam %>%
ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse() + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = OMZBase, color = "darkblue") + labs(y = "Depth (m)", x = "Particle Size Distribution Slope")
pWBPSD
PubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb < 0.5) %>% group_by(depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup()
photicBiomass <- PubDf %>% filter(depth <= 165, depth >= 155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
PubDf <- PubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBS <- PubDf %>%
ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1.2)) + geom_hline(yintercept = 160, color = "darkgreen") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = OMZBase, color = "darkblue") + labs( x = "Small particle mass (norm.)") + theme(axis.title.y = element_blank())
pWBS
#Small Particle Biomass
LubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb >= 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- LubDf %>% filter(depth <= 165, depth >=155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
LubDf <- LubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBL <- LubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 160, color = "darkgreen") + labs( x = "Large particle mass (norm.)") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = OMZBase, color = "darkblue") + theme(axis.title.y = element_blank())
pWBL
WBFig5 <- plot_grid(pWBPSD, pWBS,pWBL, nrow = 1, labels = c("B", "C", "D"))
Removed 6 row(s) containing missing values (geom_path).Removed 7 row(s) containing missing values (geom_path).
WBFig5
WBcombined <- plot_grid(WBColorMap + theme(plot.margin = unit(c(0,3,0, 3), "cm")), WBFig5, ncol = 1, labels = c("A", ""))
WBcombined
ggsave("figures/WBModelValidation.png")
Saving 8.5 x 6 in image
scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxP16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_log10(limits = c(35, 150),breaks = seq(from = 20, to = 150, by = 20)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
geom_hline(yintercept = 200, color = "darkgreen") +
theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm"))
#
#
#
# pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
# pltFlxLegend <- get_legend(pltFlx)
#
pltFlxP16
# #plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3P16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_continuous(limits = c(-1, .1), breaks = seq(from = -2, to = .75, by = 0.5)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
geom_hline(yintercept = 200, color = "darkgreen")+
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3P16
#plotly::ggplotly(plt1pos)
osms_p16 <- bds %>% filter(project == "P16") %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), group = factor(time))) + geom_point(size = 3) + geom_path() + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(-1, 1)) +
geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + labs(x = "Observed - Modeled Small Particle Flux \n µmol/m^3/day") +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = PhoticBase, color = "darkgreen")
plotly::ggplotly(osms_p16)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
#ggsave("..figures/FluxSizeShift.svg"
plot_grid(
pltFlxP16,
pltDelta3P16,
osms_p16
)
Removed 28 rows containing missing values (geom_point).Removed 28 row(s) containing missing values (geom_path).Removed 32 rows containing missing values (geom_point).Removed 32 row(s) containing missing values (geom_path).Removed 29 rows containing missing values (geom_point).Removed 29 row(s) containing missing values (geom_path).
ggsave("figures/P16FluxRelate.svg")
Saving 8 x 8 in image
ggsave("figures/P16FluxRelate.png")
Saving 8 x 8 in image
Take one profile, attenuate it, and show what it looks like
source("ModelStuff.R")
Parsed with column specification:
cols(
Cruise = [31mcol_character()[39m,
Station = [31mcol_character()[39m,
`mon/dd/yyyy` = [31mcol_character()[39m,
`hh:mm` = [34mcol_time(format = "")[39m,
`Longitude [degrees east]` = [32mcol_double()[39m,
`Latitude [degrees north]` = [32mcol_double()[39m,
`Bottom Depth [m]` = [32mcol_double()[39m,
`Pressure [db]` = [32mcol_double()[39m,
`Temperature [degrees C]` = [32mcol_double()[39m,
`Temperature 2 [degrees C]` = [32mcol_double()[39m,
`Salinity [psu]` = [32mcol_double()[39m,
`Salinity 2 [psu]` = [32mcol_double()[39m,
`Fluorescense [mg/m^3]` = [32mcol_double()[39m,
`Beam Transmission [%]` = [32mcol_double()[39m,
PAR = [32mcol_double()[39m,
`Oxygen [umol/kg]` = [32mcol_double()[39m,
`Oxygen [% saturation]` = [32mcol_double()[39m
)
scan_for_example <- bds %>% filter(project == "ETNP", depth < 500, depth > 200) %>% select(profile, depth, DFP, use_this_DFP, ospsDZ)
#loc_station = "stn_036"
loc_station = "stn_043"
loc_depth = 225
loc_prev_depth = 112.5
allDFPs <- bds %>% filter(profile == loc_station, depth >= loc_prev_depth, depth <= loc_depth) %>% summarize(DFP = prod(DFP), use_this_DFP = prod(use_this_DFP))
loc_DFP <- allDFPs %>% pull(DFP)
loc_use_DFP <- allDFPs %>% pull(use_this_DFP)
for_single_disag <- bes %>% filter(profile == loc_station, depth %in% c(loc_prev_depth, loc_depth)) %>% select(depth, lb, nnp_smooth) %>%
mutate(depth = recode(depth, `112.5` = "Shallow", `225` = "Deep")) %>% # I have no idea how to not hard code this bit
pivot_wider(names_from = depth, values_from = nnp_smooth)
with_disag <- for_single_disag %>%
mutate(Predicted_Deep = remin_smooth_shuffle(Shallow, loc_use_DFP))
#remin_smooth_shuffle(for_single_disag$Shallow,loc_use_DFP)
for_plot_disag <- with_disag %>% pivot_longer(cols = -lb) %>% #filter(lb <= 5) %>%
mutate(name = factor(name, levels = c("Shallow", "Deep", "Predicted_Deep"))) %>%
mutate(name = recode_factor(name, Shallow = "Shallow (112.5m)", Deep = "Deep (225m)", Predicted_Deep = "Predicted Deep (225m)"))
for_plot_disag %>% ggplot(aes(x = lb, y = value, shape = name)) + geom_point() + scale_x_log10() + scale_y_log10() + scale_shape_manual(values = c(1, 6, 3)) + theme(legend.title = element_blank()) + labs(x = "Particle Size (mm)", y = "Normalized Particle Abundance (#/L/mm)")
ggsave("figures/DisagExample.png")
Saving 7.29 x 4.5 in image
ggsave("figures/DisagExample.svg")
Saving 7.29 x 4.5 in image